Example of 20 hypothetical 1ha restoration plots (10,000m2, 100m*100m) distributed at random at three reefs in the Townsville region (Davies Reef, Big Broadhurst Reef, Little Broadhurst Reef). See code and “Visualise restoration plots” for workflow.
# data packages
library(tidyverse)
library(units)
library(foreach)
# spatial packages
library(sf)
# mapping packages
library(tmap)
library(leaflet)
# create function to get make polygons for plots
set_restoration_plot_centres <- function(input = NULL, width = NULL, length = NULL) {
# Calculate the coordinates of the rectangular polygon
x <- sf::st_coordinates(input)[1, 1]
y <- sf::st_coordinates(input)[1, 2]
# set parameters
x_min <- x - (width / 2)
x_max <- x + (width / 2)
y_min <- y - (length / 2)
y_max <- y + (length / 2)
polygon <- sf::st_polygon(list(rbind(c(x_min, y_min), c(x_min, y_max), c(x_max, y_max), c(x_max, y_min), c(x_min, y_min)))) |>
sf::st_sfc(crs = 20353)
return(polygon)
}
habitats <- c("Plateau", "Back Reef Slope", "Reef Slope", "Sheltered Reef Slope", "Inner Reef Flat", "Outer Reef Flat", "Reef Crest")
cairns_geomorphic <- st_read("/Users/rof011/coraldynamics/data/Moore-Arlington-20230906220745/Geomorphic-Map/geomorphic.geojson") %>%
sf::st_transform(crs = sf::st_crs("EPSG:20353")) %>%
dplyr::group_by(class) %>%
dplyr::mutate(habitat_id = paste0(
gsub(" ", "_", class), "_",
sprintf(paste0("%0", ceiling(log10(max(1:length(class)))), "d"), 1:length(class)))) %>%
sf::st_make_valid() %>%
# mutate(habitat_id=as.factor(habitat_id)) %>%
dplyr::group_by(habitat_id, class) %>%
summarise(geometry = st_union(geometry)) %>%
mutate(area = round(st_area(geometry))) %>%
filter(class %in% habitats) %>%
drop_na(class)
## Reading layer `Moore Arlington' from data source `/Users/rof011/coraldynamics/data/Moore-Arlington-20230906220745/Geomorphic-Map/geomorphic.geojson' using driver `GeoJSON'
## Simple feature collection with 3768 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 145.9182 ymin: -16.89574 xmax: 146.2582 ymax: -16.61326
## Geodetic CRS: WGS 84
# find restorable habitats
cairns_sheltered_sites <- filter(cairns_geomorphic, class %in% c("Sheltered Reef Slope", "Outer Reef Flat"))
# identify neighbours
neighbors <- st_touches(cairns_geomorphic)
# match pairs
border_rows <- purrr::map2(neighbors, cairns_geomorphic$class, ~ {
any(cairns_geomorphic$class[.x] == "Sheltered Reef Slope") & (.y == "Outer Reef Flat")
})
# Filter rows
indices <- which(unlist(border_rows))
result <- cairns_geomorphic[indices, ]
# combine outer and sheltered
cairns_restorable <- rbind(
filter(cairns_geomorphic, class %in% c("Sheltered Reef Slope")),
result
)
cairns_plots_centroids <- cairns_restorable %>%
filter(area > set_units(10000,"m^2")) %>%
filter(class %in% c("Reef Slope", "Sheltered Reef Slope", "Back Reef Slope")) %>%
drop_na(class) %>%
st_centroid() %>%
st_cast("POINT")
cairns_plot_outputs <- foreach(i=1:nrow(cairns_plots_centroids), .combine="c") %do% {
tmp <- set_restoration_plot_centres(cairns_plots_centroids$geometry[i], width=100, length=100)
tmp
}
Select layers using the layercontrol on the left, use [ ] for full-screen viewing.
set.seed(1)
cairns_plot_outputs <- cairns_plot_outputs |> st_as_sf() |> mutate(id=paste0("Plot_ID_",seq(1,n(),1)))
cairns_plot_outputs2 <- cairns_plot_outputs[sample(nrow(cairns_plot_outputs), 20), ] # subset to 20
# visualise with thematic maps
map <- tm_view() +
#tm_tiles("Esri.WorldImagery", group = "[Satellite map]", alpha = 0.5) +
#---------- all habitats--------@
tm_shape(cairns_geomorphic |> mutate(class=as.factor(class)), name = "Reef habitats", id="area") +
tm_borders(col = "black", lwd = 0.2) +
tm_fill("class", name = "Benthic habitats", Title = "Benthic habitats", id="area", palette = c("Plateau" = "lightskyblue1", "Back Reef Slope" = "darkcyan",
"Reef Slope" = "darkseagreen4", "Sheltered Reef Slope" = "darkslategrey",
"Inner Reef Flat" = "darkgoldenrod4", "Outer Reef Flat" = "darkgoldenrod2",
"Reef Crest" = "coral3"), alpha = 0.6) +
#---------- restorable habitat--------@
tm_shape(cairns_restorable |> mutate(class=as.factor(class)), id="area", name = "Restorable habitats") +
tm_borders(col = "white", lwd = 0) +
tm_fill("class", title="Restorable habitats", name = "Benthic habitats", id="area", palette = c("Sheltered Reef Slope" = "darkslategrey", "Outer Reef Flat" = "darkgoldenrod3"), alpha = 0.6) +
#---------- restoration plots--------@
tm_shape(cairns_plot_outputs2, id="plot_id", name = "Restoration plots (1ha)") +
tm_fill(col="red", alpha=0.6) +
tm_borders(col="black") +
#---------- scale--------@
tm_scale_bar(width=200)
map |> tmap::tmap_leaflet() |>
leaflet::addProviderTiles('Esri.WorldImagery', group = "Satellite map", options=leaflet::providerTileOptions(opacity=0.6, maxNativeZoom=18,maxZoom=100)) |>
# leaflet::addProviderTiles('Esri.WorldTopoMap', group = "<b> [Restoration]</b> habitats", options=leaflet::providerTileOptions(maxNativeZoom=19,maxZoom=10)) |>
leaflet::addLayersControl(position="topleft", overlayGroups=c(
"Satellite map",
"Reef habitats",
"Restorable habitats",
"Restoration plots (1ha)"),
options=leaflet::layersControlOptions(collapsed = FALSE)) |>
leaflet::hideGroup(c("Reef habitats")) |>
leaflet.extras::addFullscreenControl(position = "topleft", pseudoFullscreen = FALSE)